CONVERGENT NUMERICAL SCHEMES FOR THE 
COMPRESSIBLE HYPERELASTIC ROD WAVE EQUATION 
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Abstract. We propose a fully discretised numerical scheme for the hypere- 
lastic rod wave equation on the line. The convergence of the method is estab- 
lished. Moreover, the scheme can handle the blow-up of the derivative which 
naturally occurs for this equation. By using a time splitting integrator which 
preserves the invariants of the problem, we can also show that the scheme 
preserves the positivity of the energy density. 



1. Introduction 
We consider the compressible hyperelastic rod wave equation 

Ut - Uxxt + SuUx - j{2UxUxx + UUxxx) = 0. (1) 

The equation is obtained by Dai in [B] as a model equation for an infinitely long rod 
composed of a general compressible hyperelastic material. The author considers a 
far-field, finite length, finite amplitude approximation for a material where the first 
order dispersive terms vanish. The function u = u{t, x) represents the radial stretch 
relative to a prestressed state. The parameter 7 e R is a constant which depends on 
the material and the prestress of the rod and physical values lie between -29.4760 
and 3.4174. For materials where first order dispersive terms cannot be neglected, 
the KdV equation 

Ut + UUx + Uxxx = 

applies and only smooth solitary waves exists. In contrast, the hyperelastic rod 
equation ([T]) admits sharp crested solitary waves. 

The Cauchy problems of the hyperelastic rod wave equation on the line and 
on the circle are studied in [S] and respectively. The stability of a class of 
solitary waves for the rod equation on the line is investigated in [5] . In [T^] , Lenells 
provides a classification of all traveling waves. In [13 1 the authors establish, 
for a special class of initial data, the global existence in time of strong solutions. 
However, in the same papers, they also present conditions on the initial data for 
which the solutions blow up and, in that case, global classical solutions no longer 
exist. The way the solution blows up is known: In the case 7 > 0, there is a point 
a; € K and a blow-up time T for which liuit^T Uxit, x) = —00 for (for 7 < 0, we 
have liuit^T Ux(t, x) = 00). 

To handle the blow-up, weak solutions have to be considered but they are no 
longer unique. For smooth solutions, the energy /^(m^ + u^) dx is preserved and 
i7^(R) is a natural space for studying the solutions. After blow-up, there exist two 
consistent ways to prolong the solutions, which lead to dissipative and conservative 
solutions. In the first case, the energy which is concentrated at the blow-up point 
is dissipated while, in the second case, the same energy is restored. The global 
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existence of dissipative solution is established in In the present article, we 
consider the conservative solutions, whose global existence is established in 

There are only a few works in the literature which are concerned with numerical 
methods for the hyperelastic rod wave equation. In |13| . the authors consider a 
Galerkin approximation which preserves a discretisation of the energy. In [3], a 
Hamiltonian-preserving numerical method and a multisymplectic scheme are de- 
rived. In both works, no convergence proofs are provided and the schemes cannot 
handle the natural blow-up of the solution. 

In this paper, we propose a fully discretised numerical scheme which can compute 
the solution on any finite time interval. In particular, it can approach solutions 
which have locally unbounded derivatives (the condition Ux € L^(M) allows for an 
unbounded derivative in i°°(R)). A standard space discretisation of ([1]) cannot 
give us global solutions. To obtain these solutions, we follow the framework given 
in [TT]. With a coordinate transformation into Lagrangian coordinates, we first 
rewrite the problem as a system of ordinary differential equations in a Banach 
space (Sections [2] and [3]). We establish new decay estimates (Section!?]) which allow 
us to consider solutions defined on the whole real line. We discretise the system 
of equations in space (Section [SJ and time (Section [7]) and study the convergence 
of the numerical solution in Section (5] In Section [51 we explain how to define 
a converging sequence of initial data. This construction can be applied to any 
initial data in H^{M.). Finally, in Section [HI numerical experiments demonstrate 
the validity of our theoretical results. Moreover, the time splitting discretisation 
enables the scheme to preserve invariants and we can use this property to prove 
that the scheme preserves the positivity of a discretisation of the energy density 

+ dx, see Theorem 18.21 

The results of this paper are also valid for the generalised hyperelastic rod wave 
equation 

Ut - U^xt + ^g{u)x ~ li'^U.j^Uxx + UU^xx) = 0, u\t=Q=Uo. (2) 

However, for simplicity only the numerical discretisation of equation ([I} will be 
analysed. Equation ^ was first introduced in [3]; it defines a whole class of equa- 
tions, depending on the choice of the (locally uniformly Lipschitz) function g and 
the value of the parameter 7, which contains several well-known nonlinear disper- 
sive equations. Taking 7=1 and g{u) = 2ku + (with k > 0), equation ^ 
reduces to the Camassa-Holm equation [2]; For g{u) = 3u^, equation ^ becomes 
the hyperelastic rod wave equation ([T]); For g{u) = 2u + and for 7 = 0, equa- 
tion ^ leads to the Benjamin-Bona-Mahony (BBM) equation (or regularised long 
wave) [1]. 

2. The Semigroup of Conservative Solutions 

The purpose of this section is to recall the main results of [TT] about the con- 
servative solutions of the hyperelastic rod wave equation ([1]) . The total energy for 
the hyperelastic rod wave equation is given by the norm, which is preserved in 
time for smooth solutions. An important feature of this equation is that it allows 
for the concentration of the energy density [u^ +u^) dx on set of zero measure. To 
construct a semigroup of conservative solution, it is necessary to keep track of the 
energy when it concentrates. This justifies the introduction of the set V defined as 
follows. 

Definition 2.1. The set T) is composed of all pairs (u, such that u belongs to 
H^{M.) and is a positive finite Radon measure whose absolute continuous part, 
^ac, satisfies 

Mac = (u^ + ul) dx. 



NUMERICAL SCHEMES FOR HYPERELASTIC ROD 



3 



The measure ^ represents the energy density and the set T) allows to have 
a singular part. The solutions of ((IJ are constructed via a change of coordinates, 
from Eulerian to Lagrangian coordinates. An extra variable which account for the 
energy is necessary. Let us sketch this construction. We apply the inverse Helmholtz 
operator (Id — dxx)~^ to ([T|) and obtain the system of equations 

ut + juux + Px = (3a) 

P-Pxx='-^u' + lul. (3b) 

By using the Green function of the Helmholtz operator, we can write P in an 
explicit form, i.e., 

P{t,x)^U c-\'^-^\{l^u' + lul){t^z)dz. (4) 



2 / ju 2 2 

We also define 



Q{t,x):^Px{t,x)^~\ / sgn{x-z)e-\^-^\{^lu' + lul){t,z)dz. (5) 



2 JR 2 2 

Next, we introduce the characteristics y{t,£,) defined as the solutions of 

with 2/(0,^) given. The variable y{t,^) corresponds to the trajectory of a particle 
in the velocity field ^u. However, the Lagrangian velocity will be defined as 

U{t,0=u{t,y{t,0). 

From ([5a|l . we get 

Ut{t, = ut{t, y) + ytUx{t, y) (wf + ^uux)(t, y) = -Px{t, y) = -Q{t, y). 

As it can be checked directly from ([3]), smooth solutions satisfy the following trans- 
port equation for the energy density: 

(u2 + ul)^ + (7u(u2 + u2))^ = {u^ - 2Pu)^. (6) 

After introducing the cumulative energy H{t,^) as 

H{t,0 / {u^ + ul)dx, 



we can rewrite the transport equation ^ as 

To obtain a system of differential equations in terms only of the Lagrangian variables 
X = {y, U, H), we have to express (jH) and ([S|) in terms of these new variables. This 
can be done (see [TT] for the details) and we obtain 



Q(i,e) = -i / sgn(C-ry)c-^sn(?-r,)fe(«)-.(„))(^3_27^2y^^2i/^)(,,)d^. 



Finally, we obtain the following system of differential equations 

yt = (7a) 
Ut = -Q (7b) 
Ht^U^ - 2PU, (7c) 

which we rewrite in the compact form 

Xt = F{X). 
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The derivatives of P and Q are given by 
and 

Qdt.O^lHi+^-^U^y^-Py^. (8) 
Thus, after differentiating ([7]), we obtain 

Cct = lU^ (or y^t lU{), (9a) 



Uit = jH^ + [^^U' -P) Vi. (9b) 

H^t = -2Q Uy^ + (3C/2 - 2P) U^, (9c) 

where we denote y{£^) = C(^) + ^. 

The mapping _F is a mapping from E to E, where E' is a Banach space that we 
now define. We denote by V the space defined as 

v^ife Ct{R) I /c G L\m)}, 

where C6(R) = C(R) n L°°(R). The space is a Banach space for the norm 
ll/llv il/llL°° + ll/jllia- The Banach space E is then defined as 

E==V X xV 

with norm ||/||^ := |l/|ly + ll/ilifi + ll/lly- In IHIi the existence of short-time 
solutions of ([T]) is established by a standard contraction argument in E. The solu- 
tions of ([T]) are not in general global in time but for initial data (<^0j ^^Oi Hq) which 
belongs to the set which we now define, they are. 

Definition 2.2. The set T consists of all (C, U,H) G E such that 

{C,U,H) e \W^^°°{M.)f and lim H(0 = (10a) 

^ 0,H^ > 0, y^ + H^> c almost everywhere, for some constant c > (10b) 
y^H^ = y^U"^ + ?7| almost everywhere. (10c) 

The set F is preserved by the flow, that is, if X{Q) € F and X{t) is the solution to 
([7]) corresponding to this initial value, then X{t) & F for all time t. The properties 
of the set F can then be used to establish apriori estimates on the solutions and 
show that they exit globally in time, see [TT] for more details. We denote by St the 
semigroup of solutions in F given by the solutions of ((T]) . 

Given an initial data (it, /i) G T), we have to find the corresponding initial data 
in F; we have to define a mapping between Eulerian and Lagrangian variables. To 
do so, we set 

y(0 = sup {y I m((-oo, y)) + y < C} , (Ha) 

H[0^i-y{0. (lib) 

U{i)=uoy{i). (11c) 

We define X ~ L{u, fi) and L maps Eulerian to Lagrangian variables. When /i = /lac 
(no energy is concentrated), equation (|llap simplifies and we get 

MO 

yiO+ / {u^ + ul){x) dx = ^. 

J —oo 

Reciprocally, we define the mapping M from Lagrangian to Eulerian variables: 
Given X = {y, U, H) G F, we recover (u, /i) ~ M{X) G I? by setting 

u{x) ~ U{£^) for any ^ such that x = y{£,)i (12a) 

^i = y#{H^dO■ (12b) 
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Here, y^^iH^ dS,) denotes the push- forward of the measure d£, by the mapping y. 

In conclusion, the construction of the global conservative solutions is based on the 
change of variable from Eulerian to Lagrangian. However, this change of variable 
is not bijective. The discrepancy between the two sets of variables is due to the 
freedom of relabeling in Lagrangian coordinates. The relabeling functions can be 
identified as a group, which basically consists of the diffeomorphisms of the line 
with some additional assumptions (see Given X = (y,U,H), the element 

X o f = {yof,Uof,Hof) is called the relabeled version of X with respect to the 
relabeling function /. We can check that 

M{X) = M{Xo f), 

that is, several configurations in Lagrangian variables correspond to the same Euler- 
ian configuration. In this article, we will not be too concerned with the aspects of 
relabeling invariance. They have however to be taken into account to establish the 
semigroup property of the semigroup of solutions in Eulerian variables. We also 
use the relabeling invariance in Section IH] to construct initial data in a convenient 
way for some particular initial conditions. We define the semigroup of solutions Tt 
in Eulerian coordinates as 

Tt=AIoSto L. (13) 
Finally, we recall the following main result from |11| . 



Theorem 2.3. The mapping T: 2?xR_|_ — )- V, where V is defined by Definition \2.l[ 
defines a continuous semigroup of conservative solutions of the hyperelastic rod wave 
equation ([T]), that is, given {u,fi) £ T), if we denote by t (u(t)^fi{t)) = Tt{u^fl) 
the corresponding trajectory, then u is a weak solution of the hyperelastic rod wave 
equation 

The function y{t, ^) gives the trajectory of a particle which evolves in the velocity 
field given by ^u{t,x). If u is smooth, then it is Lipschitz in the second variable 
and the mapping ^ — > y{t, ^) remains a diffeomorphism. We denote its inverse by 
X — >■ y~^{t,x). In this case, the density p{t,x) is given by 

p(t,x) = — 3-77 TT- (14) 

We can also recover the energy density as 

{u'+ul){t,x)^^{t,y-^t,x))). (15) 

In the following sections, we design numerical schemes which preserve the positivity 
of the particle and energy densities as defined in and . 



3. Equivalent System of ODEs in a Banach space 

In this section, we reformulate the hyperelastic rod wave equation ([3]) as a system 
of ordinary differential equations in a Banach space as this was done in |11) but 
where we decouple the functions y, U and H and their derivatives y^, and 
that we denote by q, w and h. We set C(*j ■= vit^O ~ ^ ^(C; 1) = 9(*; ~ 1- 
If we assume that 



q = y^, w = and h = 



(16) 
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then (O and © rewrite 



Ct 


= yt 






Ut 


= ~Q, 




Ht 


= - 2PU, 


Vt 


= qt 


= iw, 




Wt 


2 ^ 2 




ht 


= -2QUq+(3U^ 



2P)w, 



(17a) 
(17b) 
(17c) 
(17d) 

(17e) 

(17f) 



where P and Q are given by 



P 



- ssnii^r,)iy{i)-ym {^^^U\ +lh) (77) d?7 (18) 



and 



Q 



"2 / sgn(^-77)e- 



- sgn(5-r,)(y(C)-y(,,)) ( 3-27 



2 -U'q+lh){rj)dri. 



(19) 



Note that equation (|17p is semihnear in the variables {q,w,h). Since the terms P 
and Q have similar structure, in the remaining of the paper most of the proofs will 
be established just for one of them. Now, we do not require ([T5| to hold any longer 
and, setting Y := (C, U, H, v, w, h), we obtain the system of differential equations 

Ytit) = GiY{t)), 

where G is defined by p7p . In the remaining, we will sometimes abuse the notation 
and write Y = {y,U, H,q,w,h) instead of F = {C,U, H,v,w,h). Then, we implicitly 
assume the relations y(^) = ({^) + ^ and q = v + 1. The variables y and q are the 
physical ones but do not have the proper decay/boundedness properties at infinity 
and this is why ( and v have to be introduced. The system ([T7)) is defined in the 
Banach space F, where F is given by 

F := L°°{R) X (i°°(R) n L^{R)) x L°°(R) x L^{R) x L^{R) x L^{R). 

For any Y = {(, U, H, v, w,h) ^ F we use the following norm on F: 



\Y\ 



m 



L2 



\u\\ 



\H\\ 



L2 



L2 



\h\\ 



L2 



The following proposition holds. 



Proposition 3.1. The mappings P : F ^ H^(R) and Q : F ^ H^{R) belongs to 
Ci(F,i?i(R)) andG: F F belongs to C^(F,F). Moreover, given M > 0, let 



Bm = {X eF\ ||X||^<M}. 
There exists a constant C{M) which only depends on M such that 



\\piY)\\H^ + \mY)\\H, 



-iv) 

dY^ ' 



L{F,m) 



^(F) 

dY^ ' 



< C{M) (20) 



L{F,m) 



and 



for all Y eB 



\\GiY)\\, 



-iv) 

dY^ ' 



< G{M) 



(21) 



L{F,F} 



A/- 



Here, abusing slightly the notations, we denote by the same letter P the function 
P{t,£,) and the mapping F i-> P. The same holds for Q. The norms L{F,H^( 
and L{F,F) are the operator norms. 
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Proof. First we prove that the mappings F H' P and F i—> Q as given by ([T5)l and 
(Unj) belong to C^{F,L°°{M.) D L'^{R)). We rewrite Q as 



-C{0 



2 



(22) 



where xb denotes the indicator function of a given set B. We decompose Q into 
the sum Qi + Q2, where Qi and Q2 are the operators corresponding to the two 
terms in the sum on the right-hand side of ([55]). Let /i(^) = X{4>o}(06~'^ and A 
be the map defined hy A: v t-^ v. Then, Qi can be rewritten as 



Qi 



-AoR{Ym, 



(23) 



where R is the operator from F to (R) given by 



Rivm 



C(0/3-27 



2 ' ' 2 

The mapping A is a continuous hnear mapping from L^(]R) into i^(M) ni°°(M) as, 
from Young inequahties, we have 



L2 



L2 • 



||/i*w|Il2 < \\h\\^^ \\v\\^2 and !l/i*u|li=o < \\h\ 
For any F G B^/, we have 

IIQiiL^nLo^. < C{M) \\AoR\\^,^^^ < C{M) \\R\\^._ < C{M) 



(24) 



for some constant C(M) which depends only on M. From now on, we denote 
generically by C{M) such constant even if its value may change from line to line. 
The same result holds for Q and P. Since R is composed of sums and products 
of maps, the fact that R: F ^ Li^ is follows directly from the following 
short lemma whose proof is essentially the same as the proof of the product rule 
for derivatives in R. 



Lemma 3.2. Let 1 < p < 00. If Ki G C^{F, 1°" 
the product K1K2 belongs to C^{F,LP{M.)) and 



andK2 e C^{F,LP{R)), then 



0{K,K2) (^^j^j ^ A',(y)^(y)[y] + i^.(r)^(r)[y]. 



With this lemma in hands, we thus obtain that 

dR 
dY 



dY 



and 



{Y)[Y] = e'; {^-^ (Ct/'(1 + v) + 2UU{1 + v) + U^v) + - 

< C{M). 



-iY) 
dY^ ' 



Then, Qi is in C^{F,L^ 
dQi 



ni°°(M)), 



,dR, 



Qy {Y)[Y] = _(CA(i?(r)) - A{ — {Y)[Y])) 



and 



dQ 



dY 



Hy) 



< C{M). 



L(F,L2nL«=) 
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We obtain the same result for Q2, Q and P. We differentiate Q and get 

Q,^lh+'-^U\-Pq, (25) 
see (|5]). Hence, the mapping F n- is differentiable, 

^^HY)[Y] = U+ ^-^{2UU + U'q) %{Y)[Y]q Pq, 



and 



dY ^ " ' 2 2 ' ^' dV 



dY 



(Y) 



< C{M). 

L{F,L^) 



It follows that Q belongs to C^{F,H\M.j) and w ^ ^ - C{M). The same 

result holds for P and (PU)) is proved. By using Lemma 13.21 we get that G £ 
C^iF.F) and this proves □ 

By using Proposition 13.11 and the standard contraction argument, we prove the 
existence of short-time solutions to ([T7)) : 

Theorem 3.3. For any initial values Yq = (c^Q; Uq, Hq, vq, Wq, ho) G F, there exists 
a time T, only depending on the norm of the initial values, such that the system 
of differential equations (|17p admits a unique solution in C^([0, T], -F). More- 
over, for any two solutions Yi and Y2 such that sup^gjg j^j ||li(t)||^ < M and 
suPtg[o,T] 11^2(^)11^^ < M, then 

sup \\Y,{t)~Y2mp<C{M)\\Y^{0)^Y2m\p, (26) 

te[o,T] 

where the constant C{M) depends only on M . 

Proof. The stability result is a direct application of Proposition 13 . 1 1 and Gron- 
wall's Lemma. □ 

The system of differential equations (jl7p in the Banach space F has an interesting 
geometric property: it possesses an invariant. In fact, the following quantity 

I{Y) ■.^U'^q^ + - qh 

is conserved along the exact solution of the problem as we show now. For any Y{t) 
solution of pT|) . we have 

^/(r(t)) = 2UUtq^ + 2U^qqt + 2wwt - qth - qht = ~2UQq^ + 2U^qjw 

+ 2w(^^h + {^^^U^ - P)q'j - iwh - q{-2QUq + (SJ/^ - 2P)) = 0. 

(27) 

Additionally, we have 

Lemma 3.4. The following properties are preserved (independently one of each 
other) by the governing equations (jl7p 

(i) q, w, h belongs to L°°{M.). 

(ii) qh = C/2q2 + j^yj ^ 

(iii) qh = U'^q'^ + w'^ (or I{Y) = 0) and q > 0, h > 0, q + h > c almost 
everywhere for some constant c > 0. 

(iv) The functions y, U and H are differentiable and y^ = q, ^ w and 

= h. 

Proof. We consider the short time solution given by Theorem 13.31 on the time in- 
terval [0,T]. Given an initial data Fq which satisfies (i), then we use Gronwall's 
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Lemma and the semi-hnearity of (|17dp - (|17fp with respect to q, w, h and obtain 
that 

ikWIli- + \Ht)\\L^ + \\Ht)\\L^ < ciWqoh^ + !|wolL=e + ||/iolL=e) 

for a constant C which only depends on j|lo!li?- have already seen that I{Y{t)) 
is preserved, so that the condition (ii) is satisfied. We consider a fixed ^. Let us 
denote T = sup{i G [0, T] \ [q + h){t, ^ > for aU te[0,t]}. Since 

qh^q^U^ + w^ (28) 

for all time t E [0,r], the product qh is positive and therefore g > and /i > for 
t E [0, T"]. From the governing equations pT|) . we get 

_d ^ „ \w\+q 

dt^q + h' - + 

for some constant C which depends only on the norm of the initial data. Since q and 

h are positive, it follows from that |u;| < ^{q + h) and therefore — < 

C 1 1 ~ — 
. Gronwall's inequality yields (t) < e^^ . Hence, we have T = T 

q + h q + h 90 + ^0 

and 

C-^^{qo + ho) < {q + h){t). 
In particular it implies that T ~ T and there exists a constant c > such that 
q + h > c ioT almost every ^ and t E [0, T]. Thus we have proved that Y{t) satisfies 
the condition (iii) . The last property follows from |TT] , where it is proved that there 
exists a unique solution to the system 

Ut = -Q, (29) 
Ht = U^ - 2PU 

in the Banach space V x iJ^K) x V, where V = {f E L°°{R) \ E L'^{R)}. By 
differentiation of we obtain that {y,U, H,y^,U^, H^){t) satisfies P7)) . There- 
fore by uniqueness of the solutions, for an initial data Yq which satisfies (iii), we 
obtain that the property (iv) is satisfied. □ 

Having a closer look at Lemma [231 we now define the following set. 

Definition 3.5. The set Q consists of the elements {y,U, H, q,w, h) E F which 
satisfy the conditions (i), (iii) and (iv). 

As a consequence of Lemma 13. 4[ the set Q is preserved by the system. For any 
initial data in Q, the solution of p7)) coincide with the solutions that are obtained 
in In particular, we prove in the same way as in that 

Tiieorem 3.6. For initial data in Q , the solutions to (jl7p are global in time. 

We denote by St the semigroup of solutions to ([TT]) in Q. Thus, we slightly abuse 
the notations, as St was already introduced in the previous section, see (|13p. but, 
as we explained, the two semigroups are essentially the same. Note that global 
existence can only be established for initial data in Q and do not hold in general 
for initial data in F . 

4. Decay at infinity 

The terms P and Q, as given by (fT5)) and (fT91) which appear in the governing 
equations (|17p are global in the sense that they are not compactly supported even 
if Y is. Consequently the set of compactly supported functions is not preserved 
by the system. However, we identify in this section decay properties which are 
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preserved by the system. We denote by F"^, the subspace of F of functions with 
exponential decay defined as 

F" = {Y e F \ q,w,h e L°°{R), e^^^U, el«lw e i^(M), e^^^h e L\R)}. 

We define tlie following norm on F'^ 



\Y\ 



WW 



\h\\ 



II liirrll II lil I 

+ ||C 2 U\\l2 + ||C 2 w\\l2 + 



Given a > 1, we denote by F", the subspace of F of functions with polynomial 
decay defined as 



(i + |C|)t[/, (i + |f|)t^eL2(]R), 

{i + \^\rh eL\R)}. 



F" ^{Y eF \ q,w,he L°°{R), 

We define the following norm on F" 

\\Y\\p^=\\Y\\j, + \\q\\^^+\\w\\^^+\\h\\^^ 

+ + im^uh. + 11(1 + ICD^HU^ + + \^\rhhr. 

Theorem 4.1. The spaces F*^ and _F" are preserved by the flow of (jl7p . Consid- 
ering the short-time solutions given by Theorem YS.'Sl we have that 

(i) IfYo e F^-, then sup,g[o,T] \\Y{t, < C, 

(ii) IfYo e F'^, then sup^^^^.T] IWit, ■)\\p^ < C, 

for a constant C which only depends on T and \\Yf)\\pc (case (i)) or T and Hlolli^a 
(case (ii))- 

Proof. Let us prove the case (i). First, we establish bounds on the solutions. 
By applying the Cauchy-Schwartz inequality, we get 









2 


e 


2 



which implies that Uq G L^(R) and ||C/o|Ili < C for some constant C which depends 
only on ||el^lC/Q Similarly we get that wq € L-'^(M) and llwolj^i < C for some 

constant C which depends only on ||el^lz/;Q(^)||^i. We denote generically by C such 
a constant, which depends only on T and HyollF"- From Theorem 13.31 and Lemma 
13.41 we get that 



\Qit,-)\\ 



\witr)\\ 



\hit,-)\\ 



< C. 



ll^oa -r ll'^V'-J VIIL^ ^ ll'H'-J JUL" 

By following the same argument as in the proof of Proposition 13.11 from to 

((24)1 . but, instead, using the Young inequality || 

that 

||Q(f,.)ILi<c(||Mt,-)ILi + i) (30) 

for a constant C which depends only on ||y(t)||^e and, therefore, only on |jlollj;-= 
and T. The same estimate holds for P, that is, 



\\Pit,-)h.<Ci\\hit,-)h.+l). 

Let us denote 

J{t) := \\Uit,-)\\^, + \\w{t,-)\\^, + \\h{t,-)h.. 
From the governing equations (fT7|) , after using and (PTjl , we get 

Jit) < J{0)+C + C f J{T)dT. 



(31) 
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Hence, by applying Gronwall's Lemma, we get that, for t G [0, T], 
J{t)^\\U{t,-)\\^, + \\w{tr)h. + \\h(t,-)\\^, <C 
for another constant C. Let L{t) denotes 

L{t) = ||cl-'IC/2(t, Ol^^ + ||cl«lu;2(t, .)||^^ + ||cl«l/i(t, •) 
From the definition of Q, we get that 

Q{t, 0<C I c-l«-"l + h){t, r,) drj 
Jr 

so that 

el«lg(t,e)<C / el«le-l«-''le-l''lel''l(t/2 + /i)(i,^)d^ 
< CL{t) 

because '•^^'^ therefore 

|el«lQ(i,-)||^^ < CL{t). 

Similarly, we get that 

|el«lp(t,-)||^^ < CL{t). 
From the governing equations pT|) . we get that 



(32) 
(33) 

(34) 



(35) 
(36) 



el«lf/^(t,0 



< 



< 



< 



2cI«IQ[/(t,-) 



dT 



c'-Q{t,-) ^^||[/(r,.)ILi dr 



+ C / /(r)dr, 



(37) 



by using the apriori estimates ([5^ and ([55]) . From (fT7)) . we also obtain that 



el«l/i(t,C) 



< 



Cl'l/IQ 

t 



LI 



+ / (2 cl«lQ(r,.) ^ ||[/(r,.)ILOrfr+ 



(C el«l[/2(r,.) , + cl«lp(r,.) |1 w(r, .)ILO 







which, after using the estimates p2p . ([55)1 and yields 



< 



el«l/io 



+ C + C I{T)dT. 



(38) 



Similarly we get that 



el«lu;2(i,e) 



< 



cl«l/i(r,0 



< 



+ / (C el«IC/2(r,-) +C 
el«lw2 +C + C f I{T)dT. 



dT 



cl«lp(T,-) ||w(T,-)||ii)dr. 



(39) 



After summing dST]), dMl) and ([391), we get L{t) < L{0) + C + C L{t) dr and 
the result follows by applying Gronwall's inequality. We now turn to case (ii). We 
introduce the quantity 

m = !!(! + + 11(1 + nrw-'it, + iki + icd-m*, oiil- 



12 



COHEN AND RAYNAUD 



From (|M|1 . we get 

{l + \^\rQ<C I {l + \C\re-\^-'H^ + \v\r"{l + \ri\nu'q + h)dr^. (40) 



Since < \^-v\ + \v\ < (l + IC - ?7l)(l + b7l), we have (1+|C|) < 2(1 + |^ - 77|)(1 + |,7|) 
and 

(1 + i^ir <2"(i + i^-r;ir(i + Hr. (4i) 

Then, it follows from (gOl) that 

(l + |eirg<C / c-\i-^\il + \^-r^\ril + \^\nu\ + h)drj 

< C \\e-''{l + |z|)"||io= K{t) < CK{t) (42) 
so that 11(1 + |C|)"Q|1l=o < CK{t). We have to estimate ||(1 + \£,\TQ\\li. We have 

I1(1 + ICI)"QIIli < / (l + |^|)"e-l«-"l(l + H)-"(l+H)"(C/29 + /i)dr;dC 



(1 + |77 + z|)"e-l^l(l + |?7|)-"(1 + |77l)"(C^''Z + h) dTjdz 
<2" [ (l + |z|)"e-l-'l(l + |7/|)"(C/2g + /i)drydz (by gl])) 
< C / (1 + |z|)"e-l^l dz [{1 + IvinU^ + h) dij 



< CK{t). 

Hence, 

11(1 + iei)"QILinL== < Ci^(t) 
and the same bound holds for P. From the governing equations, we obtain 

||(l + l^l)"t/'(i,0||M < ||(l + l?irt^o||,i+ ri|2(l + iei)"0C/(r,.)IL. dr 



(43) 
(44) 



<||(l + |^|)"f/o1|,.+2 / \\{l + nrQ\r,-)\\^,dT 



+ 2 / 11(1 + |ei)"C/^(r, OIL, dr 



< ||(l + l?irt/ol,,+C / A»dr^ 



by (UH), as IIQjl^oo < C, see In a similar way, one proves that 



1(1 + \i\rw\t,Oh. < 11(1 + iei)"^^o|Li + C I K[r) dr 



and 



li(l + |e|)"Mi,e)llLi < ||(l + |ei)"/iolLi +C r/v(r)dT 
so that ^ 

A'(t) < A"(0) + C I K{r) dr 
Jo 

and the result follows from Gronwall's Lemma. 

For later use, we note that, in this proof, we have established that 



+ 



<C{\\Y\\^.) 



and 



11(1 + iei)"Qii™i + 11(1 + iei)"^ii™. < c(i|y|!^„) 



□ 

(45) 
(46) 
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for some given increasing function C, see ([55)) . (|36p and 

5. Semi-Discretisation in space 

The first step towards a discretisation of ([TT]) is to consider step-functions. We 
consider an equally-spaced grid on the real line defined by the points 

where A<^ is the grid step and z = 0, ±1, ±2, . . .. We introduce the space 

= {Y e F : each component of Y consists of 

piecewise constant functions in each intervals [CiiCi+i)}- 

The system ([T7|) does not preserve the set Fa^ of piecewise constant function. Thus, 
we define 

oo 

PAdYm-^ E piYm)x[iu,+Mo, (47) 



QAdYm^ E Q(^)(eoxK.,«.+,)(^) (48) 

i— — oo 

and consider a second system of differential equations 

Ut = -QAi 
Ht = U^- 2Pa^U 
qt = jw 

ht = -2QAiUq + {3U^ - 2PAi) 



(49) 



w, 



or, shortly, 

Yt{t)^GAdY{t)). 

Like in the preceding section, we show that this system of differential equations 
possesses a short-time solution, an invariant and that it solution converges to the 
solution of p7)) as — 0. In the next theorem we prove, by a contraction 
argument, the short-time existence of solutions to (02]). 

Theorem 5.1. For any initial value Yq = {yo,Uo, Hq, q(),WQ, ho) G F, there exists 
a time T, only depending on \\Yq\\p, such that the system of differential equations 
(|49p admits a unique solution in C^([0,T],F). 

This theorem is a consequence of point (i) in the following lemma. 

Lemma 5.2. The following statements hold 

(i) The mapping Ga^ : F ^ F belongs to C^{F,F) and 

dGAi 



dY 



-{Y) 



< C{M), (50) 

L(F,F) 



\\GAdY)h 

for any Y e Bm ■ 
(ii) For any Y Cz F, we have 

\\G{Y)~GAdY)h<C^^ (51) 
for some constant C which only depends on 
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Proof. For any function / G iJ^M), let P(/) be the function defined as P(/)(^) = 

oo 

/(6)XK,,Ci+i)(C)- Thus, we can rewrite QA^iY) and Pa^{Y) as 

i——oo 

QAdY) = P[Q(n] and PAdy) = P[P{Y)]- 

Let us prove that P is a continuous mapping from H^{R) to i°°(R) n By 
using the Sobolev embedding theorem of H^{M.) into L°°(R), we get 

l|P(/)IL» < <c||/i|^, 

for some constant C, so that P is continuous from H^{M.) into L°°(M). The 
norm of P(/) is given by 

oo 

i— — oo 

We have, for all ^ € that 

/te)' = /(C)'-2^V(r?)/«(r])dr; 

which, after integration over [Ci,Ci+i)! yields 
Hence, 

iip(/)iii. <(i + Aoii/i!i. + Aeii/diI. 

and the mapping P is continuous from H^{R) to L^(M). Since Qa^ and Pa? are 
compositions of a continuous linear map P and a map, they are also and 

^(?) = p(f(v)|fi) 

for all y € P. The same holds for Q so that (|50p follows from Lemma [3.21 Let us 
prove point (ii). First we note that ([5T|) follows directly from the definitions of G, 
Ga? and the estimate 



\\Q{Y) - QAaniL.nLoo + ll^(^) - PAdY)\\L^nL^ < CVAC- (52) 

Let us prove (|52p . We estimate ||Id — P||^(^i ^.^ni^)' '^here the norm here is the 
operator norm from H'^{R) to i°°(M) n ^^(IR). Let us consider / e H^iR), we have 



11/ - P(/)IL=. < sup 11/(0 -/(e.)ll 



For any ^ e [e.,6+i), we have |/(0-/(6:)l < %/Se||/dli2, by the Cauchy^ 
Schwartz inequality. Hence, 



||/-P(./)IL.. < v/A^ll/dL. < VAeil/llHi 
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We have 



1 1 ? ^ 



< fl{r,)dv r^\^-i.)di 



fi{v)dr,. 



Hence, 



11/ -P(/)||^.<^ 11/11^. (53) 

and we have proved that ||Id — P\\j^2f^i^oa < Cy/A£^ for some constant C. Then, we 
have 



\\Q{Y) - QAdY)h.nL- < C^A^ \\Q{Y)\\j,, < C ^ A^ 

for another constant C" which depends only on One proves in the same way 

the same estimate for P and thus we obtain (1521). □ 



Concerning our new system of equations (j49p . it is not difficult to show in the 
same way as in (P7)) that 

/Ac(r) ■.= U^q^ + ^ qh 

is also a conserved quantity along the exact solution of our problem. The system 
(Uni) is introduced because it allows for a space discretisation of the original system 
([T7]). Indeed, the set of piecewise constant functions is preserved: 

Lemma 5.3. The set F/^^ is preserved, that is, ifYo G andY{t) is the solution 
of (gni) with initial data Yq, then Y{t) £ Fa^ for all t £ [0,r]. 

The proof of this lemma is straightforward. We can now compare solutions of 
(|49| and of the original system (jl7p . 



Theorem 5.4. Given Af > and Fo,^o,A{ G F. Let Y{t) be the short-time 
solution of p7p with initial data Yq and YA^{t) be the short-time solution of (|49p 
with initial data io,A5 in the interval [0, T] . If we have 

\\Y{t)\\p < M and ||YAc(i)llF ^ f'^^ * ^ fO'^l' 
then we also have 

\\Yit) - YAdt)\\p < {\\Ya - ro.Adl + CrVA|)c^^ for all t € [0,T] (54) 
with some constant C which depends only on M . 

Proof. The proof of this theorem is a consequence of Lemma l5.2l and of Gronwall's 
Lemma. We have 

Yit) - YAdt) = Yo- yo.A? + / (G(r(T)) - GAdYAdr))) dr 



= Yo- Yo,Ai + / (G(r(r)) - G{YAdr)) + G{YAdr)) - Ga4(1a5(t))) dr 
Ja 

which yields, after using Proposition 13.11 and Lemma 15.21 

\\Yit) - YAmip < 11^0 - Yo^AiWp + C f \\Y{t) - YAdr)\\p dr + GT^, 
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for some constant C which depends only on M . Then, ([M)) fohows from Gronwall's 
Lemma. □ 



Lemma [231 and Theorem 14. 1 1 show that there exist properties of the initial data 
that are preserved by the system ([TT)). The same results - with the exception of 
property (iv) in Lemma [3.41 - hold for the system This is the content of the 

following theorem. 

Theorem 5.5. We consider an initial data Yq ^ F and the corresponding short 
time solution Y{t) of given by Theorem \5.l[ 

(i) Ifqo, Wo, ho belongs to L°°(R) then 

sup (lk(t, OIIl- + \Ht, Ollioe + \\h{t, ■)\\^^) < c 

te[o,T] 

for some constant C which depends only on T and HyollF"- 

(ii) // we have qh = U^q^ + for t = Q (or I{Yq) = Oj then this holds for all 

[0,T]. 

(iii) // we have qh = U'^q^ + w"^ (or I{Y) = 0) and q > 0, /i>0, q + h > c 
almost everywhere for some constant c > 0, then the same relations holds 
for all t e [0, T] . 



(iv) IfYo e F^, then 



ifYo e F", then 



sup \\Y{t,-)\\p. 

tG[0,T] 



< a 



sup \\Y{t,-)\\p^<C, 
te[o,T] 



(55) 
(56) 



where the constant C depends only on T and j|yo|li?e; and T and ||yo||^„, 
respectively. 

Proof. The system (|49p is obtained from pT)) by simply replacing P and Q by Pa^ 
and Qaj as defined in (|T7|) and (|^5|) . Therefore, the proofs of points (i), (ii) and 
(iii) in Lemma 13.41 which do not require any special properties of P and Qj apply 
directly to . After introspection of the proof of Theorem 14.11 we can see that 
in order to prove (|55p. we need to prove that the estimates (pO| . (|3T|) . (j35p . ((36| . 
which hold for P and Q, also hold for Pa^ and Qa^, namely, 



and 



PAdtr) <CL{t), 



(57) 
(58) 



where L{t) is defined in ([55)) and C is a constant which depends only on T and 
li^olli^c- We denote generically by C such constant. In the same way that we 
obtained ([55)) . we now get that, for any / S VF^'^(R), 









P(/)(«)l di = 




d£. 
















\fdv)\ dv 



and therefore 



i/-p(/)ii^, <Aeii/d 



(59) 
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We obtain, after using successively ([SH)) . (P5|) and ([?T|) . that 

IIQAdlii <IIQA«-QILi + ||giLi 

<Ae||Qdl^, +1) 



2 



-U q — Fq 



LI 



C(||/.|| 



LI 



1) 



<C||P||^,+C7(||/i||^,+l) 



<C(||/i|U, + l). 

We handle in the same way ||Pa^||^i and this concludes the proof of ((57)) . For any 
^ € M; we have ^ S for some i. Then, 



_ „l«|-|«iLIC.I/ 



by ([55)1 and, therefore, ||el^l(5A^(i, OH^ae !i CL{t). Similarly, we obtain the corre- 
sponding result for Pa? so that (|55)) is proved. Again, after introspection of the 
proof of Theorem 14. 11 we can check that, in order to prove ([55)) . we need to prove 
that 



11(1 + |CI)"QAc(t, ■)h^nL^ + 11(1 + l^irPAiit, ■)h^nL^ < CK{t) 
We have 

ii(i + icirgAc(t,-)iLoo < \\{i + \£,\rQ{t,-)\\^^<cK{t) 

by dH]). Since e^"'' < e'^«e«'"'' for any (^,77) e we get 

(i + ieirQA«(t,-)iL. < E 



(60) 



(l + |C|)"c-l«--"l(C/'9 + /i) dTjd^ 

(i + ieire-i«-''i([/2g + /i)d?7de 



<e-« E 

= c^^ ( [ (l + |C|)"e-l«-"l(C/2g + /i)d?7d^<CA'(i), 

by (|^5)) . The corresponding results for P are established in the same way and this 
concludes the poof of (|60)) . □ 



In order to complete the discretisation in space, we have to consider a finite 
subspace of Pa^- Given any integer N, we denote R = A^A^ and we introduce the 
subset Pr of P defined as 

FR = {YeF : 

UiO = QiO = HO = HO = 0, for ah ^ G (-c», -P) U [P, ex.), 

C(0 = Cco, P(e) =i^oo, foraUee [P,oo), 

C(0 = C-00, P(0 = for all e e (-c», -P), 

where C±oo and Pqo are constants}. 

The set P^ basically corresponds to functions with compact support ([/, q, w and 
h vanish outside a compact set). We do not require that the functions C and H 
have compact support (C and H belongs to L°° with no extra decay condition) but 
we impose that they are constant outside the compact interval [— P, P]. We denote 
P{A^,R} = Pr n Pa^. The set P{A{,fl} is not preserved by the flow of (j^^ because, 
as mentioned earlier, P and Q do not preserve compactly supported functions. 
That is why we introduce the cut-off versions of P and Q given by 

P{Ai,R}{Ym= E ^(n(6)XK.c.+,)(0, 
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(61) 



JV-l 
i=-N 

and define a third system of differential equations 

qt = 7""^, 

ht = -2Q{A?.i?}C/g + (3[/' - 2P{Ae,fl})«^, 

or, shortly, 

It is clear from the definition that the system (pT|) preserves -F{A4,_r} and therefore, 
since -P{A4,_r} is of finite dimension, the system (|5T|) is a space discretisation of p7)) 
which allows for numerical computations. To emphasize that we are now working 
in finite dimension, we denote 

Ct = C{Ai,R}{t,^i), Ui = ?7{Ac.fl}(^,^t) and so on for Hi,qi,Wi,hi, and for 
i = {-N, . . . , - 1}. We have 



-N 



Again, we can show that 



I{Ai.B}(Y) ■■^Ufq^+W^-q.h, (62) 

are conserved quantities along the exact solution of problem Finally, note 

that -F{Aj.i?} is contained in F"^ and F". Concerning the exact solution of we 
have the following theorem. 

Theorem 5.6. For an initial values Yq ~ (yo, C^Oi ^ZOi w;0: ^o) G F? there exists 
a time T, only depending on the norm of the initial values, such that the system of 
differential equations (j6ip admits a unique solution in C^{[0,T], F). 

This theorem is a consequence of point (i) in the following lemma. 

Lemma 5.7. The following statements holds 

(i) The mapping G{A4,_r} ■ F F belongs to C^{F,F) and 



\G{AiMY)\ 



dY 



-{Y) 



< C{M), (63) 

L{F,F} 



for any Y £ Bm ■ 

(ii) For any Y £ F"^ , we have 

\\GiAi,R}{Y) - GAdY)\\p < Cc~", (64) 

for some constant C which only depends on \\Y\\p,,. 

(iii) For any Y € F"' , we have 

\\G{Ai,R}iY) ~ GAdY)\\p < C(VAe+ (65) 



for some constant C which only depends on \\Y\\pa 
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Note that for Y{t) solution of (pT|) . we have 

sup \\Y(t,-)\\p, <C and sup \\Y{t, ■)\\p.^ < C, 
te[o,T] te[o,T] 

where C depends on ||Fo||^e and ||yo||p„; respectively. This follows from (|55p . 
dMl) and dMl)- 

Proof of Lemma^ For any function / e i°°(R) n ^^(M), let Pfl(/) be the func- 
tion defined as Pi?X/)(0 = f{£,)X[-R,R)- Thus, we can rewrite Q{A(,,r}(X) and 
P{Ai,R}{Y) as 

Q{AtR}iy) = PRlQAdY)] and P{Ai,R}iy) = Pfl[PA?(F)]. 
The operator P/j is a projection from L°°(K) fl L^(R) into itself and therefore its 
norm is smaller than one. Hence, ([55)) follows from ([50)) . Let us prove (ii). We 
consider Y G F'^. We have to prove 

\\Q{AiM}iY) - QAcCnlL^nio. + \\P{Ai,R}{y) - PAdY)L.r^L^ < Ce-«. (66) 
By (gSj), we have ||cI«IQ||^^ + ||cl«lp||^„ < C. Hence, 

\\Q{Ai.R.} - QacLoo = sup |Q(^,)I < C sup c-l«'l = C7e-«. 

l€i|>-R l«il>-R 

We have 

|C,|>_R |jAe|>_R 

and therefore ||Q{Ac.fl} — '9a5||^2 ^ Ce~^. We prove in the same way the corre- 
sponding result for P and it concludes the proof of ([55)) . The estimate ([M)) follows 
from ([55)) . Let us prove (iii). We consider Y £ F". We have to prove that 

\\Q{AtR}iY)-QAdY)\\^.^^^ 

+ \\P{M,R.}iY) - PAdY)h.r.L^ < C(v/A|+ -^). (67) 
By we have ||(1 + leD^QILoo + 11(1 + |^|)"^1Il=c, < C. Hence, 

\\Q{Ai,R}^QA4L^^ sup \Qm<C sup (l + |e.|)-" = C(l + i?)-". (68) 

We have 

||Q{AC,i?,} - QAc||i2(R) = ll<3A?|li2(K\[_fl,fl]) 

< IIQa? - Q\\l^{m\[-r,r]) + II'3IIl2(r\[-_r,_r]) 
<C(v/Ae+||0|L2(R\[_^,^])), (69) 

from ([52]). Since 



(1 + ii?i)"" / (i + mrQ'dc 

Jm\\-R.R] 



'm\[-R,R] 
<C(l + i?)-", by (HI]), 

the estimate ([63 follows from ([68l) and ([Ml)- □ 



Again, the system ()6ip preserves properties of the initial data: 

Theorem 5.8. We consider an initial data Yq Cz F and the corresponding short 
time solutionY{t) of (|6ip given by r/teorem 15.61 Then, Y{t) satisfy points (i)-(iv) 
as given in Theorem \5.5\ 

Finally, for any initial data in Fo G F"^, resp. Yq e F", we obtain the following 
error estimate for bounded solutions. 
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Theorem 5.9. Given Yq and Yq a^.b. i'n F'^, let Y{t) and y{A^,R}{t) be the short- 
time solutions of (jl7p and (j6ip . respectively, with initial data Yq and Yo^a^,r, 
respectively. If we have 

\\Y{t)\\p, < M and ||F{AC,i?}(t)||^ < M for all t € [0,r], 
then we have 

sup ||r(i,-)-nAC,i?j(i,-)L < C^(ll^o-ro.Ac,7?||^ + v/Ae + c-^), (70) 

tG[0,T] 

where the constant C depends only on AI . For Yq and Yq r in F°' , we have that 
\\Yit)\\p^ <M and ||y{A«,i?} < for all t e [0,T], 

then 

sup ||y(t,.)-r{A5,fl}(t,-)L<c'(||yo->o,A5,fl||^. + VAe + ^). (71) 

te[o,T] ^ ' 

Proof. We have 

\\Y{t, •) - F{A?,i?}(i, OIU ^ 11^0 - YoA^-nWp + 

' \\G{Y{T,-))^G{AiMY{AtR}(.rr))\\pdT. (72) 

By Proposition 13.11 and Lemmas 15.21 and 15 . 71 we get 

||G(y(r,.))-G{AC,fl}(nA«.i?,}(r,-))ll^^ 

<||G(r(T,.))-GA«(F(r,-))||^ 

+ ||GAe(r(r,-))-G{A«,i?}(r(r,.))||^ 

+ ||G{A«,fl}(n^,-))-G{AC,ii}(i^{A5,i?,}(T,-))||^ 

1 

< G((A02 +e-^+ ||r(r,.) -r{A«,ii}(r,-)||^) 

for a constant G which depends only on M. Hence, ([7D|) follows from ([7^ after 
applying Gronwall's Lemma. The proof of ([7T|) is similar. □ 

6. Approximation of the initial data and Convergence of the 
Semi-Discrete solutions 

6.1. Approximation of the initial data. The construction of the initial data 
Yq^A(,r is done in two steps. First, we change variable from Eulerian to Lagrangian, 
that is, we compute Yq £ Q such that X = {yQ,Uo, Hq) £ satisfies 

Uo^UQoyQ. (73) 

In the new set of variables, we can solve ([TT)) or, rather, its discretisation ((6T|) . Note 
that, given uq G iJ^(M), there exists several Yq G Q such that ([75]) holds. This is 
a consequence of relabeling invariance and this fact will be used in the numerical 
examples of Section [U] Here, we present a framework valid for general initial data 
in iJ^(M). In Section [51 we define the mapping L from V to For uq e i7^(R) 
and fiQ absolutely continuous, it simplifies and reads 

rvaiO 

2/o(0+ / ("0 + < J ^, (74a) 
UQ = UQoy and i^o=Id-?/o- (74b) 

Then, we set 

90 = ^0,5, w = l7o,5, h = Ho^^. (74c) 
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As earlier, we denote vq = I — qa and Co = Id — yo- We have 

^090 = 90^0 + 'W'o' Qo + ho = 1, Qo > 0, ho > for almost every ^ e M. 

(75) 

The element Yq = {yo,Uo, Ho,qo,wo, ho) belongs to Q. The second step consists 
of computing an approximation of Yo in i*{Aj,i?} • In the following theorem, we 
show how the change of variable given by ((71)) deal with the decay conditions. For 
simplicity, we drop the subscript zero in the notation. Let us introduce the Banach 
spaces H^''^ and iJ^'" as the subspaces of with respective norms 



and 







2 




2 




e 2 u 


+ 


e 2 -u^ 








L2 


L2 



1(1 + lel)^«f.. + 11(1 + 



Theorem 6.1. Given u and Y as given by (|74p . we have 

(i) u e H^'" if and only ifY eF", 

(ii) ue H^^'^ if and only ifY e F''. 

Proof. Let us assume that u E H^'^ . By definition, we have h = {v? + u^) o yy^. 
Hence, 



S^hiOd^^ I c\^Hu' + ul)oy{OydOd^ 
\y'''-^^\{u' +ul){x)dx 
\v''(^)~^\c\''\{u^ +ul){x)dx 



< ell!'(?)-«lli,- / c^'^^u^ + ul){x)dx <oo. 



Using (|75)) . we get 



Jm Jm 
In order to prove that J^e^^^U^ d£^ is finite, we decompose the integral as follows: 



{«6R|9<^} 



el«IC/2d^ 



{CeR|g>i} 



el«IC/2 d^. 



We have 



{CeKk<^} 



el«IC/2d^ < \\U\ 



{«eK|<?<^} 



el«l d^ 



<\m] 



{CeR|/i>i} 



e'^' d^, a,s q + h 



and 



< 2 ||[/||^^ / /iel«l di<C j e^^^hdi < oo 

J{^es.\h>^} Js. 

[ e\^\U^di<2[ 
J{£eR|«>i} J{. 



< 2 

< oo 



{«eR|g>i} 1 

e^^^qhdS,, as < qh bv ([75l). 

UeR|g>i} 
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Hence, J^e^^^U^ < oo. Let us now assume that Y G F"^. Then, 



P^^^\h{£_)d^ 

< / el^(«)-«lel«l/i(Od^ 
Jr 



and uq E H^ '^. The case (ii) is proved in the same way. □ 

As a consequence of this theorem and Theorem 14. 1[ we obtain 

Theorem 6.2. The spaces H^ '^ and H^ '^ are preserved by the hyperelastic rod 
equation: If uq G H^''^, then w(t, •) G H^-'^ for all positive time and, similarly, if 
Uq G H^'" , then u(t, ■) G i?^'" for all positive time. 

To the best of our knowledge, these decay results are new, even for the Camassa- 
Holm equation (case 7 = 1). They have to be compared with where it is 
established that the only solution which has compact support for all positive time 
is the zero solution, i.e.. the compactness of the support (which is a kind of decay 
condition) is not preserved by the equation. 

Let us now construct the approximating sequence for the initial data. From (j75p , 
we get that 

0<g<l, 0<h<l 

and 

U^ = w<^<^{h + q) = ^. (76) 



Given an integer n, we consider and R such that - = + 



so 



that ri — > 00 if and only if A^ — > and i? — > 00. We introduce the mapping 
Ia^ : — > which approximates functions by piecewise constant functions, 
that is, given f G L'^, let 

and set 

00 

i— — 00 

We define F„ = (y„, f/„, i?„, (7„, z«„, /i„) as follows. Let 

As usual, we denote q = 1 + v and qn = 1 + Vn- Let us define the weighted integrals 

It' i^ud^ 



We set 

Un{0= E for* = -A^,...,iV-L 

i=-N 

We define 



HniO^plJ hndv) if ^e[-R,R] 
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and Hn{^) = /i„ drj ii ^ e (-00, -R), i7„(C) = K dyy if ^ G (i?, 00). For ?/„, 
we set 

and y„(0 = ^ " irn(--R) if C G ("Oo, -i?), y„(C) = ^ - i?n(i?) if ^ e (i?,oo). The 
definition of P is given in the proof of Lemma 15.21 The following theorem states 
that Yn approximates Y in -FjA^.j?} scud satisfies additional properties which will be 
useful in Theorem l8.21 where we prove that the positivity of the energy is preserved 
by the numerical scheme. 

Theorem 6.3. Given Y ^ Q , there exist a sequence Yn € -F{A5,fl} such that 

lim ||r„-y|l^ = 0, (77a) 

n— f 00 

and 

qnhn > U^Qn + Qn + hn — 1, for all 71 > and for all ^. (77b) 

Moreover, we have 

\\Yn\\p, < and \\Yjp^ <C\\Y\\p„ (77c) 

for Y F'^ , resp. Y G F" , and where the constant C which does not depend on Y 
and n. 

Proof. Let us first prove (|77bp . Since q + h = 1 (see (j75p ). we obtain g„ + h„ ~ 1 
from the definitions of Vn (recall that g„ = 1 + «„) and /i„. We consider a fix given 
interval / = [^i, ^i+i] and, for convenience, denote by an integral without boundary 
the weighted integral / f{^)dS, ~ //'^^ /(C) so that, for ^ ^ I, Qn ~ J qd^, 
Wn = J w d^ and hn = J hd^. Using Jensen's inequality, we get that 



ql + Ulql + wl = ( j qd^y + Ul{ J q d^Y + ( J w d^Y 
<Jq' d( + U^I q'd^ + J w^dS, 
= Jq'd( + Ul Jl'd^ + I (^(l - 9) - q'U') 



^qn + Kjq'd^-J iq'U') d^. (78) 
Using the Cauchy-Schwarz inequality and the definition of Un , we obtain 

Hence, yields 

<ll + ulql + ^1 < 9" 
which, as qn + hn = 1, is equivalent to qnhn > Unql + Wn- Let us now prove (j77ap . 
A direct computation shows that 

||PfllA?(/)llL2 < ll/IL., (79) 

for any / G L^(K) and any n. Since lim„_j.oo ||PfllA5(/) — f\\p2 = for any smooth 
function / with compact support, we obtain, by density and (j79p . that the same 
result holds for any / e L^(M). Hence, 

lim \\qn — q\\p2 = 0, lim \\wn — w\\p2~0 and lim ||/i„ — /i||^2 = 0. 



On the interval / ~ we have 

JqHv){u{v)~u{0)dv 



J q^ drj 



- 2 
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as \U^\ < \, see ([75|) . Hence, ||[/„ — t^|licx,(_^ ^) < ^ and 

< + I|C^IIl-((-oo -_R)U(fl.oo)) • (80) 

Since [/ e i/^M), lim^^ioo C^n = and (gOj) yields lini„_,oo \\Un - U\\^^ = 0. We 
have 

^ E TFTT^ / / q'{v){U{0-U{v))'d^dv, (81) 
after applying Cauchy-Schwarz. For ^, ?y e /, we have 



{UiO - U{f^)Y = ( / U^{fi) dfjy < Ae / U^ifj)' drl < AC / t/| d^. 
Hence, (gT]) yields 

iic/„-p;^x^)iii. < E ^i^^r r^^'^^ /^"V|df;<(Ao^iit/^ii^,. 

It follows that 

\\Un - c/|L. < ||c/„ - PR{u)h. + \\u - Pr{u)\\^. 

<A^\\U^\\^, + \\U\\ 

L2((~oo,-_R)U(/?,oo)) 

and therefore lini„_j.oo \\Un — U\\]^2 = 0. The function h belongs to L^(M) because 
h ~ h'^ + U^q"^ +10^, by (|75|) . A direct computation shows that 

llPfllAc(/)IL. < ll/IL., (82) 

for any / G i^(M) and any n. Since lini„^.oo ||P_rIa^(/) — /IIli = for any smooth 
function / with compact support, we obtain, by density and (|82p . that the same 
result holds for any / € L^(R). Hence, lim„_j.oo \\hn — ^H^i = and therefore 

lim \\Hn-H\\^^ =0. 

Since yn = £. — Hn and y = ^ — iJ, we get also that lim„_>.oo Wun — vWl^ = 0. Let 
us look at the bounds on the decay of Y . We assume Y € F'^ . We have 



1 ^ /-C.+i /-f.+i 



= E /"^' / "'el«le-l"lel"l|Mr?)Nr;de 
< c'^^ ^ / el"! |/i(?7)| dr? < 3||yj| 
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after assuming, without loss of generality, that < 1. Similarly one proves that 
/g e'^'w^ < C||F||pe- It remains to estimate J^^U^e^^^ d^. For any 77,^ S 
we have 



Hence, 



2 I !t'em{ri)dr^ 



i.n 



< — — T (by Cauchy-Schwarz) 

- J^:^W{ri)dr, 

<u\0 + I^^^^\u' + w')i0d^ 



for any e Then 

^-1 f6+ 



< E r"(ei«i(c^'(o+ r\u'+w'md^)d^ 



< E [^'^\u'+w'me\^\did^ 

< (l + 2A^c^«)||y||^. <7n|^.. 

Thus we have proved that |jl^n||j;^e < C|jy|jj;,Q for a constant C which does not 
depend on Y and n. One proves in the same way that < C ||y|j^„. □ 

6.2. Convergence of the Semi-Discrete solutions. Let Y{t) and 1{AC,/S}(i) be 
respectively the solution of p7)) with initial data Yq and the solution of (pT|) with 
initial data lo.Af.fl. We assume Yq £ F^. Given T > 0, we consider the fixed time 
interval [0,T]. Since Yq g Q, the solution Y{t) exists globally and 

sup \\Y{t,-)\\p. <M 
te[o,T] 

for a constant M which depends only on T and ||yo||j;,c, see Theorems 13.61 and 14.11 
The solution Y{A5,fl} does not necessarily exist globally in time. However, we claim 
that there exists n > such that for any and R such that A^ + < we have 

sup \\Y{ACR}it,-)\\ < 2M. (83) 
te[o,T] 

It implies in particular that the solution Y^^^n^ is defined on [0,T]. Let us as- 
sume the opposite. Then, there exists a sequence A^^, and < T such that 
limfc^oo A^fc = 0, limfe^^oo i?fe = 00, 

sup ||r{A?,i?j(i,-)|| - 2M. 

tG[0,tfc] 
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From ([70)) . we get 

sup ||r(t, •) - y{Ae„fl.}(t, ■)\\p< C{M){\\Yo - Yo.ac.mJIf + TA? + e-^*"). 

(84) 

The constant C{M) depends on M but not on A^fe and Rk- Thus, we have 

2Af= sup < rfe,-)ll + ||>^(ifer)-nA?..H.}(^fc,-)|| 

te[o,tfc] 



which leads to a contradiction as the right-hand side in the last inequality above 
tends to M when k tends to infinity. Once is established, Theorem 16.41 follows 
from (|70[). The same estimates can be obtained for Fo G i^". Without loss of 
generality, we assume that the approximating sequence satisfies \\Yq — io,A5,i?||^ ^ 

'"^M^ where C{M) is given in ([84)1 . so that FjA^.A} exists on [0,T]. Then, we have 
the following theorem. 

Theorem 6.4. Given Yq G F'"' , for any T > 0, there exists a constant n > such 
that, for all and R such that A^ + < , w^e have 

sup ||r(t, •) - y{AC,i?}(t, Ollf < C( llFo - i"o,Ae,i?||j. + + e"^). 

tG[0,T] 

The constant C depends only on Iji^oll^^e O'l^d T. Correspondingly, given Yq g F°' , 
we have 

sup \\Y{t, ■) - Oil;. < C( llFo - yo.ACflllj. + ^) 

tG[0,T] -t ' 

and C depends only on \\Yi)\\pc, and T. 

7. Discretisation in time 

In this section, we deal with the numerical integration in time of the system of 
differential equations (|61l) which corresponds to the semi-discretisation in space of 
system (jl7p . The flow of this system of differential equations has some geometric 
properties and it is of interest to derive numerical schemes that preserve these prop- 
erties. Such integrators are called geometric numerical schemes, see for example the 
monograph [9]. Thus we will look for numerical schemes preserving the invariants 
([S^ of our system of differential equations. Moreover, this last property will enable 
us to show that the numerical schemes preserve the positivity of the energy den- 
sity. These invariants are quartic functions of Y and we are not aware of schemes 
preserving quartic polynomials, this is why we first split the system of equations 
(|6ip into two pieces. Each sub-system will then have quadratic invariants and we 
can use a numerical scheme preserving these invariants. The following sub-systems 
read 



Ci,t 


= 






= 






= 




li.t 


= IWi 








C-^'^Uf P.)g.. 


hi,t 


= (3f/f 





i = -N,...,N~l (85) 



or shortly 



Yt = G,(Y), 
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where Y{t) = (K[A^,_R}(i7 a-^d similarly for Gi. We also define the system 
of differential equations 

Cm = 7t/« 
Hu = Uf - 2P,U, 

g,,t = i = -N,...,N-l (86) 

Wi,t = 

hi,t = ~2QiUiqi, 

or shortly 

i^t = G2(y). 

The space -^{ACii} is finite dimensional. We denote F = R^^^^. The mapping 
from F to F{ac.h} 

is a bijection, where we define 

N-l 

CiO = (CjX[c.,c,+i)(0) +C-ivX(-oo,-i?](0 + CwX[_R„oo](0 

and similar definitions for the other components of Y. This mapping is in addition 
an isometry if we consider the norm 



l^lli2(K2") ll^lu° 



H 



l"lli2(R2N) + ||w||;2(^2N) + 1 1 ^ 1 1 j2 (]a2N ) 7 (87) 



where 



(AC E 



for any z £ R^^. In the remaining, we will always consider the norm given by 
([87)) for F so that the bounds found in the previous sections directly apply. In 
particular, we have the following lemma, which is a consequence of Proposition 13. II 
and the same arguments that lead to Lemmas 15.21 and 15.71 

Lemma 7.1. The mappings Gi : F F and G2 ■ F ^ F belong to C^{F, F) and 

dGi 



and 



\\G2{Y)\\p + 
for any Y G Bm, where 



BY 
dG2 



dY 



< C{M), 



L{F,F) 



< C{M), 



L(F,F) 



Bm = {Y(^f\ ||y||^<M}. 

As this was done in the last sections, one can show that both systems posses 



Qih,, see dHll), as first integrals. That is mY)Gk{Y) = 



for all y, for fc = I, 2 and for i — —TV, . . . , — 1. In particular, this implies that 
every solutions of (|HS|) or dHS]) satisfy IiiY{t)) ^ h{Y{0)) for i ^ -N, . . . ,iV - 1 
and < > 0. Having a closer look at the differential equations ([55)1 and one 
sees that the invariants are now quadratic functions ([7 is constant for (|85p and q is 
constant for ) and we therefore use a numerical scheme that preserves quadratic 
invariants. 
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Proposition 7.2. Let us apply a Runge-Kutta scheme with coejficients satisfying 

biAij + bjaji ~ bibj for all i, 7 = 1, . . . , s (88) 

to the system (j85p . then it conserves exactly the invariants Ii{Y) = Ufqf+w^ — qihi 
for i = —N, . . . , N — 1. The same holds if we apply the scheme to (|86p . 

Proof. The proof of this proposition is a simple adaptation of the proof of Theorem 
2.2 from [HI Chapter IV]. Let us start with system ([55)) . Dropping the indexes and 
the bars for ease of notations, we first write the invariant I{Y) as 

I{Y) ^ Y'^D{Y)Y + d{YfY 

/O \ 





C/2 -1/2 

1 

yo -1/2 / 

For the Runge-Kutta method, we write Yi —Y^ + h bjKj with Ki = 6*1(^0 + 

hJ2j=i o-ijKj). From the definition of the method, of the matrix D{Y) and of the 
vector d{Y), it follows that 

s s 

I{Yi) ^ Y^D{Yi)Y, + d{YifYi = {Y^ + hY^b.Kif D{Yo){Yo + hY,bjK,) 



with Y = (C, U, H, q, w, h), D{Y) = 



and d{Y) = {F. 



i=l 



Y;[D{Yo)Y^ + hY^b,KfD{Yo)Y„ + hY^b,Yo^D{Yo)K3 



i=l 



+ h^Y. brb,KlD{Ya)K,. 
Writing Ki = Gi{Yi) with Yi = Yq + hJ2'j=i ^ij^j^ obtain that 

s 

I{Yi) = YjD{Yo)Y„ + 2hY,b^Y^D{Yo)Gl{Y,) 

s 

+ h^Y. (^'^J - ^'"y ~ bjaj,)KjD{Yn)Kj. 

ij = l 

The last term in the above equation vanishes due to condition (|88l) . By definition 
of the problem and of the matrix D{Y), we have D{Yq) = D{Yi) because U is 
preserved and since I{Y) is a first integral for we get Y-^ D{Yi)Gi{Yi) — 0. It 
thus follows 

/(Yi) - Y;[D{Yo)Y^ + - /(Fo) 
and the Runge-Kutta scheme applied to conserves the invariant liY). 

/O 0\ 
q2 



The proof for system ([BS]) is similar, take D{Y) = 



d{Y) = (0,0,0,0,0,-g)^ 







1 

yo 0/ 



and 



□ 



Let us consider the following differential equation yt{t) = f{y{t)). The implicit 
midpoint rule 

y, = y, + Atf{yi±y^) 
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satisfies tlie condition (|88p and thus preserves quadratic invariants. The imphcit 
midpoint rule will be the building block for the construction of the schemes we 
will use for the numerical experiments in Section [21 For other schemes preserving 
quadratic invariants, we refer to [9] for example. 

As a direct consequence of Proposition 17. 2[ we have the following result. 

Theorem 7.3. Let us apply a Runge-Kutta scheme <&Xt; fesp. with coejfi- 

cients satisfying (|88p to the system (|85p . resp. (j86p . with time step size At. Then 
the Lie- Trotter splitting 

has order of convergence one and preserves all the invariants Ii for i = — iV, . . . , TV— 
1. The Strang splitting 

*At $Xt/2 ° *At ° *it/2 

is symmetric, has thus order of convergence two and preserves all the invariants It 
fori^ ~N,...,N - 1. 

If we take for i = 1,2, the implicit midpoint rule, we obtain a first order 
splitting scheme for (|6ip that preserve exactly the invariants (a second order scheme 
is obtained using the Strang splitting). This will be the schemes that we will 
consider in the numerical experiments of Section [21 

8. Full discretisation 

Our concern is now to combine the results from the last two sections and to 
show that our numerical schemes converge to the exact solution of the system of 
equations (jl7p . We integrate Y(t) on the time interval [0, T] and obtain Yj for the 
time steps jAt,j = 0,..., Nt where At = We have the following convergence 
result. 

Theorem 8.1. Given initial values Yq in F'^ and Yq S Fh, for the Lie-Trotter 
splitting we have 

1 5, At (Fo ) - At (?o ) 1 1 ^ < C ( 1 1 Fo - I'o 1 1 p + + e-^ + At) , (89) 



max 
je{o,...,NT} 

where we recall that St stands for the semigroup of solutions to pzp and, where the 
constant C depends only on HYolli^e, ||^o||j^c andT. Correspondingly, given initial 
values Yq in F°' and Yq G Fr, we have 

max \\SMY^)-'^3U%)\\p<C{\\YQ-Y4p + ^i+-^+At), (90) 

where the constant C depends only on \\Yo\\pa, ||Yo||^„ and T. The same results 
hold for the Strang splitting with second order accuracy in time, that is, when we 
replace At with At^ in (US]) . 

Let us denote Y{t) = St{Yi^) and 

^ , ((j + l)At ~ mM%) + - jAt)<&(,-+i)At(^o) 

for t e [jAt, (j + I) At]. We can rewrite as 

max \\St{YQ) - '^t{%)\\p < C{\\Yo - Yo\\p + yA^ + e"^ + At). 

Proof of Theorem 18.11 To estimate the total error 

\\SjAt{Yo)^'i>jAt{Yo)\\p 

we split it in time and in space. Let us start with the error in time. The proof 
follows basically the steps of the standard proof of the convergence of numerical 
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scheme for ordinary differential equations. Tlie crucial point is that we guarantee 
here that the convergence rate in time is independent of the discretisation step in 
space. Let us first prove the following claim: Given M > 0, for any Y G Bm and 
Z € Bm, we have 

+ At (Gi(F) - GiiZ) + G2{Y) - G2iZ)) + 0{Af), (91) 

where (pAt{Z) stands for the exact flow of (pi]) at time At with starting values Z. 
Here, and in the following, the O-notation stands for an element in F satisfying 

\\0{e)\\p<CiM)e 

for all e > 0, where the constant C(M) depends on M but is independent on R 
and on the space grid size A^. We first show that the midpoint rule 



$i,(r) = F + AtG, 



applied to equation resp. is at least first order accurate. To do this, let 

us introduce the mapping K : F x F F given by 

K{Z,Y) = Z-Y-AtG^{^±^). 
We have K{^\^{Y),Y) = 0. Since 

dK - Atacwz + y 

^(r) _ < C{M) (by Lemma [771]), there exist C(M) such that, for At < 



and 

dK 



-jjfj-f^, we have that ^-{Y) is invertible. By the implicit function Theorem, we get 
that ^'Xt(^) is well-defined. Moreover, also following from the implicit function 
Theorem, we get that 

\\<i>iAY)\\p<CiM). 

Then, 

*L(?)^v- + A«5.(? + f<;, (£Lffl±I 

^Y + AtGi{Y) + 0{At^) 
by Lemma [7. II Using Lemma [77T] again, we obtain for the exact flow of (|55|) that 

^it{Z) = Z + AtGi(Z) + 0{At^). 
Following the same arguments, we obtain that 

*L(*At(n) = ^AtiY) + AtG^i^itiY)) + OiAt^) 
and for the composition of the exact flows 

^AtifAtiZ)) = vltiZ) + AtG^ivlm) + 0{Ae). 
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Hence, 

^ Y - Z + AtiGi{Y) - GiiZ)) 
+ At (G2(r + AtGiiY) + 0{Ae)) - G2{Z + AtGi{Z) + ©(At^))) 
+ 0{At^) 

^Y~Z + At{Gi{Y) - Gi{Z) + G2{Y) ~ G^iZ)) + OiAt"). (92) 
We consider now the splitting error. We have 

^At(^) -Z^ AtG{Z) + 0{At^) 



and 

and thus 
Hence, 



ip\t{Z) -Z = AtGi{Z) + 0{Ae) 
vltivltiZ)) = v\t{Z) + AtG^ivltiZ)) + 0(A<2). 



vltiAtiZ)) - ipAtiZ) = AtG{Z) - AtG,{Z) 

- AtG2{Z + AtGi{Z) + 0{Af)) + 0{Ai^) 
= At{G{Z) - Gi{Z) - G2{Z)) + 0(At2) 
= 0{Af), (93) 

a.s G = Gi + G2- Combining and (|92p , we obtain (pij) and the claim is proved. 
Let us now set 



M= sup ||^t(yo)||^. 

te[o,T] 



For a given At, we define 



(94) 



jAt = max{j e {0, . . . , A^T - 1} I ||$jAt(^'o)|U < 2M for all j < j}. 
For j < 7 At, we get from that 

\\'^U+i)At{Yo) ~ nj+i)At{Yo)\\p < {l+C{M)At) ||$o)At(?o) - ^u)At{Yo)\\p+0{At^). 
By induction, it follows that 

i 

\\^U+i)At{Yo) - 'Pu+iMYo)\\p < \\OiAt^)\\ ^(1 + C(A/)Ai)' 



fe=0 



< \\0{At^ 



1 



C(A/)At 



and therefore 



*0-+i)At(^'o) = ^(j+i)At{Yo) + OiAt). 



(95) 



We claim that there exists a constant C{M) such that for all At < (j^^, we have 
jAt = Nt — 1 and therefore ([M]) holds for all j < iVr — 1. Let us assume the 
opposite. Then, there exists At^ such that limfe^.oo Atk = and jAtk < — 1- 



By definition ([M)) . we have 



(jAU+l)Atfc 



_ >2M. Then, dHS]) implies 



2M < 



^OAt,+l)AtJ>o) - '/5(j^t,+l)AtJ>o) 



VU&t^^+l)Atk{Yo) 



< 0{Atk) + M 
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which leads to a contradiction when k tends to oo. Finally, for the total error in 
space and time, we have: 

\\SjAtiYo) - ^jAti%)\\p < \\SjAtiYo) - ^jAtin)\\p + - <^>MYo)\\p , 

where all the functions are evaluated at time jAt for j < Nt- The first term can 
be estimate using Theorem l6.4l and we thus obtain 

^^^max^ ^ \\S,At{Yo) - <PjAti%)\\p < C{\\Yo - Yo\\j, + ^ + e"^). 

For the second one we use and this concludes the proof of the theorem for the 
Lie- Trotter splitting. If we had taken the Strang splitting instead, we would have 
obtained an error in time of order two since this scheme is symmetric. The proof 
for initial data in F" is the same. □ 



Our next task will be to show that our schemes preserve the positivity of the 
particle density and of the energy density as does the exact solution of (fT7|) with 
initial data given by Theorem 16.31 In order to prove this result, we introduce F°° 
defined as 

= {Y={y,U,H,q,w,h) e F \ \\q\\^.^ + \\w\\^^ + \\h\\^.^ < 00} 
with the norm 

\\Y\\F.^^\\Yy + \\q\\^^ + \\w\\^^ + \\h\\^^. 

We know that the space F°° is preserved by the governing equations p7)) . see 
Lemma [231 Using the semilinear structure of (jl7dp - (jl7fp with respect to q, w, h, 
one can show in the same way that ([?!]) was shown, that, for a given M > 0, 

< C{M) (96) 

L(F°°,_F°=) 

for any Y e Bj^^ ^ {Y £ F°° \ \\Y\\p^ < M}. The same resuh holds for the 
mappings Ga^, Ga(,,r, Gi and G2. In particular we can prove, as in Theorem 18. II 
for the proof of (I9T|), that 

'^At{Y)-^At{Z)=Y - Z 

+ At (Gi(F) - Gi{Z) + G2(r) - G2{Z)) + 0{At^), 

where the definition of C'(-) is replaced by 

\\0{e)\\p^<G{M)e. 

Here, F°° — F ~ K^^^^ but equipped with the norm derived from jHl^cx,, see (|57|. 

Theorem 8.2. We consider an initial data which satisfy 

q^h'^>{U°q^r + {w^)^ q°>0, h'^ > and q^ + h'^ > c 

for all i = —N, . . . , N —1, for some constant c > 0. Then, given T > 0, there exists 
n> 0, which depends only one, \\Y^\\prx, cmdT, such that if, A£^ + -^ + At < ^,the 
positivity of the particle density 1/q and of the energy density h are preserved by 
our numerical discretisation, that is, 

ql > and h^ > 0, 

for i = -N, ...,N-1 and j = 1,...,Nt- 



\G{Y)\\ 



-iY) 
dY^ ' 
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Proof. The main idea of the proof is to control the growth of l/{q^ + h^). To 
do so we adapt the proof of Lemma 13.41 to this discrete situation. Let M = 
2 sup ||vt(^o)||j?3e- As in the proof of Theorem 18.11 we can prove that for At 

tG[0,T] 

small enough (the bound depending only on M), we have 

\\'^kAti%)\\p^<2M 
for all fc = 0, . . . , Nt- For k < Nt, we have, by definition of our scheme, that 
1 1 a,f+^ - a? + f^+' - 



1. 



k+l 



k+i 



(3(C/f )2 - 2P{Y''))wf) + OiAt^) 



Hence, using the bounds we get 
1 1 



< 



AtC{M) 


( 




+ 


9f 


+ At 













q- ' - + h^' qf + K 
Let us prove by induction that, for At small enough (depending only M), 



1 



1 



(97) 



(98) 



for i = —N, . . . , N — I, all k = 0, . . . , Nt and where C(M) is the constant given 
in (P7|) . By definition of our initial data, these assumptions hold for fc = 0. We 
assume now that (j98p holds for fc = 0, . . . , j and we want to prove that it also holds 
for J ' + 1. We set M = 1.^^c:{M)t ^ ^ Since the numerical schemes preserve the 
invariant qfh^ = {U^q^Y + (^f )^: obtain in particular that 



for all fc = 0, . . . , Nt- From this, it follows that 



(99) 

h^) as > and 



h'^ > 0. For k < j, we get from (j97p and our induction hypothesis that 



1 



1 



q- ' - + h'i+^ qf + h 
From the above equation, we get 
1 



< 



AtC{M) 



.^k+l 



,k+l\ 



1 



A/At I 



(100) 



k+i 



k+l 



< 



and therefore 



< 



< 



We have 



lim 



1 - 2C(Af )At - MC{M)At^ 



1 

(1 - 2C(A/)At - MC{M)At'^y 
1 

c(l - 2C(A/)At - A/C(M)At2)sT ' 

_ 1 2C(A/)T 



< M. 



At^o c{l-2C{M)At- MC(M)At^)±t c 

Therefore, by taking At small enough, depending only on the value of M and not 
on the number of induction steps j, we get 



+1 



< M. 



34 



COHEN AND RAYNAUD 



Using the above inequality and pOOp . we obtain 

1 1 - . V 1 



g^ + z^r- g^ 

so that j^i^^j^i > for a sufficiently small At. By we have that qi^^hj^^ > 
and therefore 

> and h{+^ > 0, 

which concludes our proof by induction. □ 

Now we go back to the original set of coordinates. Given an initial data uq G 
H^'^iR) or iJi'"(M), we construct the initial data Yq as given by (EH). Then the 
function u{t, x) defined as 

u{t,x) ^U{t,^) for y{t,^)=x (101) 

is well-defined, is a weak solution to ([3|) which corresponds to the global conservative 
solution. The definition (|10ip of u{t, x) means that for any given time t the set of 
points 

(2/(t,0,t/(i,0) eK' for^eK 

is the graph of u{t,x). Let i = + + At so that n tends to infinity if and 
only if A^, At tend to zero and R tends to infinity. We consider an approximating 
sequence Fo,n which satisfies the conditions (j77ap and (j77cD of the sequence of initial 
values which is constructed in Section |6l Let Yn{t) ~ $(yo.„). From Theorem 18.11 
we obtain the following convergence theorem. 

Theorem 8.3. The full discretised scheme provide us with points which converge 
to the graph of the exact conservative solution u(t,x). Indeed, if uq G H^''^{M.), we 
have 

max \iynitJ,^^),Unitj,C^)) ^ {y{tj,Cd,U{tj,^,))\ 

1— — N,...,N—l 
3=0, ...,Nt 

<c(||yo-^b||^ + yAe + c-^ + At), 

where the constant C depends only on |lMo||^i,e and, if uq G H^'"{M.), 
max \{yn{tj,C,),Unitj,^i)) - iyitj,^i),U{tj,^,))\ 

2— — A/,...,A/ — 1 
j=0,...,NT 

<C{\\Yo~Yo\\^ + ^C+j^ + At), (102) 
where the constant C depends only on ||mo||^i,c,- 



' <MA<C(M)(l + -L + MA<) 



Since 

|y(i,em)-y(i,6)l 



q{t,£.)di 



<CAe, 



where C depends only on ||yo|li?cxi, we have an apriori upper bound on the density 
of points of the graph of u we can approximate by our scheme. 

In the case where uq does not belong to i7^'"(M), we can approximate mq 
by functions Mo,fc G 7J^'"(R), which converge to uq in i/^(R). From [Tl], we 
know that the change of variable ([7^ produces sequences Fq.A; and Ip such that 
limfc_j.o ||i^o,fe ~ ^11^ = 0. In this way. by using the results done for functions in 

, we can approximate the exact solution Y{t) and prove convergence. However, 
since Hlo.fcllp-o is not uniformly bounded with respect to k, we lose the control on 
the error rate (the term -^hri) which is given by (|102p . 
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9. Numerical experiments 

In this section, we present some numerical experiments for the hyperelastic rod 
wave equation ([T]). In order to demonstrate the efRciency of our schemes, we will 
numerically compute three types of traveling waves with decay, see Figure [T] The 




Figure 1. Traveling waves with decay with speed c = 1: smooth 
(7 ~ 0.2), peakoii (7 = 1), cuspon (7 = 5). 

derivation of the cusped (7 > 1), resp. smooth (7 < 1), solutions follows the lines 
of [12]. We refer for example to [Hj for a thorough discussion on the peakon case 
(i.e. 7 = 1). 

9.1. Smooth traveling waves with decay (7 < 1). According to the classifica- 
tion presented in |12) . for a fixed 7 7^ 0, traveling waves u{x — ct) are parametrised 
by three parameters, M, m and the speed c. Moreover, they are solutions of the 
following differential equation 

(^-^ (103) 

c — -yu 

For positive values of 7, a smooth traveling wave with decay with m = inf^jgu u{x) 
and M — ma'Xxes.^ix) is obtained if z ~ m < M < c/7, where z := c — M ~ m. 
For our purpose, we have to set to = so that the solution decays at infinity. This 
gives us the conditions c = M and 7 < 1. We thereby obtain the initial values 
for our system of differential equations (pT|) by solving (jl03p numerically. To do 
this, some care has to be taken as u i~> ^ F{u) is not Lipschitz. We instead solve 
Uxx = F'{u)/2. Once this is appropriately done we get the initial values [/q = u, 
Wo = Ux- We then set yo = 90 = 1, '^o = Uq + and = J^"^ Iiq. These 
initial values do not correspond to the ones defined by (j74p but they are equivalent 
via relabeling and one can check that (yo, Uq, Hq, qq, wq, Hq) e Q. Figure [5] displays 
the exact solution together with the numerical solutions given by the ODE45 solver 
from Matlab, the explicit Euler scheme, the Lie- Trotter and the Strang splitting 
schemes at time T — 7. We plot the points 

{y{tA^),Uit,C^)), for*=-iV,...,7V-l, 

which approximate the graph of the exact solution u(t, x) for t ^ T. The initial 
value is a smooth traveling wave with parameters 7 = 0.2, to = 0, A/ = c = 1, 
see Figure [TJ We took relatively large discretisation parameters = 0.25 and 
At = 0.1. We observe that the explicit Euler scheme gives a less accurate solution 
than the other schemes. We also observe that, even for these large discretisation 
parameters, the splitting schemes have the same high as the exact solution, thus 
following it at the same speed. We do not observe any dissipation. Since both 
splitting schemes give relative similar results, in what follows, we will only display 
the results given by the Strang splitting scheme. We finally note that all schemes 
preserve the positivity of the particle density but only the splitting schemes conserve 
exactly the invariants from Section [7] (these results are not displayed) . We finally 
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Figure 2. Exact and numerical solutions of a smooth travel- 
ing wave with decay. Solid line=exact, dashdotted line=ODE45, 
Stars=Explicit Euler, Square=Lie- Trotter, Diamond=Strang. 



want to mention that for negative values of 7, smooth traveling waves with decay 
also exist. They are obtained if c/7 < m = Af < z. 

9.2. Peakon (7 = 1). The Camassa-Holm equation, i.e. equation ((TJ with 7 = 1, 
possesses solutions with a particular shape: the peakons. A single peakon is a 
traveling wave which is given by 

u(i,.T) =ce~l^^'=*l. 

We note, that at the peak, the derivative of this particular solution is discontinuous. 
The initial values are then 

go = 1, ho^UQ + wl, Ho ^ ho{ij) drj 

J —00 

In Figure [31 we display the numerical solutions given by the explicit Euler scheme 
and the Strang splitting for a single peakon traveling from left to right with speed 
c ~ 1, see Figure [TJ For readability reason, we do not display the solution given by 
the ODE45 solver, but we note that this numerical solution is very similar to the 
one given by the splitting scheme. Due to the discontinuity of the derivative, we 
have to take smaller (in space) discretisation parameters: = 0.05 and At = 0.2. 
We note more grid-points before the peak and very few just after it, but the speed 
of the wave is still relatively close to the exact one. As in the preceding case, only 
the splitting schemes preserve exactly the invariants of our problem. 

9.3. Cusped traveling waves with decay (7 > 1). Let us now turn our atten- 
tion to cusped traveling waves. For 7 > 0, according to the classification given in 
[12| . cusped solutions with 0/7 = nisxxi£mu{x) and m = ini^^j^u^x) are obtained 
if z = 7Ti = < c/7 < A/. This gives us the condition c = M and thus 7 > 1. The 
cuspon u(x) satisfies (jl03p . which yields for the indicated values of the parameters 

u^ = -^Fiuj^-(^^^)\ (104) 

for a; > and with the boundary value at zero given by u(0) = ^. For such 
boundary value, the differential equation (jl04p is not well-posed and the slope at 
the top of the cuspon (that is a; = 0) is indeed equal to infinity. However, we can 
find a triplet X = (y, U, H) in F which corresponds to this curve, that is, such 
that -I- u^dx) = M{X), see ((T^) for the definition of the map M. Due to 

the freedom or relabeling, the representation of the curve (x, u{x)) is not unique: 
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-5 



10 



Figure 3. Exact and numerical solutions at time T = 5. Solid 
line=exact, Stars=Explicit Euler, Diamond=Strang. 



For any difFeoniorphism (<p(^), u{ip{^))), we obtain an other parameterization of the 
same curve. Here, we look for a smooth (^(^) (and we set = (p{£,)) such that 
U = u{(p{^)) = u{y{^)) is smooth, even if u is not. We introduce the function 

, , r dz 



Since 



dx 
du 



/F{u) 



by p04p . if we choose 



u{0 = - - e, yiO = giUiO) 

7 

then we get, at least for ^ e [0, a triplet for which U{£,) — u{y{£,)). We set the 
energy density by using (jlOcp and get 



However, in this case, 



c — 



M -U U 



1 



so that (0) = oo and it is incompatible with the requirement that all the deriva- 
tives in Lagrangian coordinates are bounded in L°°(R), see (|10ap . Thus, we take 

In this case, we have 



2/7 



re 
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and 

i/aO)-^(Af7-c)* 
T 

is finite. The problem we face now is that the functions are given only on the 
interval [0, ^) and lim^_j.ji y{£) = oo. We know that the tail of the cuspon behaves 

as u{x) ~ ^e~V?=^ as x tends to oo, see [12]. Since we require that — ^ 

remains bounded, we would like to have U{C) sa :^e~v^? for large ^. Therefore we 
introduce the following partitions functions xi ^-iid ^2 defined as 

f 1 if e < a 

xm = \-T^A^~b) for [a, 6] 
[o ifx>& 

and X2(0 = 1 — Xi, where a <h are two parameters. We finally set 

7 

and 

By a proper choice of the parameters a and 6, we can guarantee that y^{£,) > for 
all ^ > 0. We extend X{£^) — (y{^), U{£,), H{£^)) on the whole axis by parity and we 
obtain an element in T such that (jl2p is satisfied. Figure 0] displays and U{(,). 
Figure [5] displays the exact solution together with the numerical solutions given by 

0.5p 

0.4 

0.3 

0.2 

0.1 
* 
-0.1 
-0.2^ 
-0.3 
-0.4. 

Figure 4. The function y{^) (left picture) and the function U{^). 
Note that these functions are smooth while uq{x) is not Lipschitz, 
see Figure [TJ 



7 

H^ = U^y^ + ^. 




the explicit Euler scheme and the Strang splitting scheme at time T = 6. As before, 
we note that the numerical solution given by the ODE45 solver is very similar to 
the one given by our splitting scheme. The initial value is a cusped traveling wave 
with parameters 7 = 5,m = 0,M = c = 1, see Figure [TJ For the discretisation 
parameters, we take = 0.1 and At = 0.1. We see that, even for initial data with 
infinite derivative Ux{0) = ±00, the spatial discretisation converges. For the time 
discretisation, as expected, explicit Euler is less accurate than the other schemes. 
We also remark that only the splitting schemes preserve the positivity of the particle 
density and conserve the invariants. We finally note that, for negative values of 7, 
an anticusped traveling wave with c/7 = Taiiix^B_u{x) and m — supx^g^u{x) is 
obtained if c/7 < m = A# < z. 
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Figure 5. Exact and numerical solutions of a cusped traveling 
wave with decay. Solid line=exact, Stars=Explicit Euler, Dia- 
mond=Strang. 

9.4. Peakon-antipeakon collisions. In Figure[5]we display a collision between a 
peakon and an antipeakon for 7 = 1. For this problem, the initial value is given by 



The numerical solutions are computed with grid parameters = 0.1 and At = 0.1 
until time T = 8. Once again we notice that the spatial discretisation converges. 
Let us now see what happens for a peakon-antipeakon collision with 7 7^ 1. In 
Figure [7] we present a similar experiment as the above one, but where we use 7 = 6 
and T ~ 2. Here, we plot the graph given by the points 



for t — T. From the right part of Figure [7] we see that only the splitting schemes 
preserve the positivity of the energy density. As always, only the splitting schemes 
conserve exactly the invariants. 

9.5. Collision of smooth traveling waves. We want now to study the behaviour 
of the numerical schemes when dealing with a collision of smooth traveling waves, 
as this in an important feature of our numerical scheme to be able to handle such 
configuration. To do so, we consider the following initial value 



Figure [5] displays the exact solution (i.e. the numerical solution with very small 
discretisation parameters) for 7 = 0.8. It is remarkable to see that even for such 
solution, our scheme performs very well. In order to get a better understanding 
of this problem, we look at the evolution of the waves with time. Figure [H] shows 
this evolution together with a zoom close to the collision time. We now present 
the results given by the numerical schemes with grid parameters = 0.25 and 
At = 0.1 in Figure [TOl We have also checked that only the splitting schemes 



u(0,2;) = c-l^l -c^l^^il. 



fori = -iV,...,iV-l 



u(0, x) 
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X 

Figure 6. Peakon-antipeakon collision for 7 = 1. Stars=Explicit 
Euler, Diamond=Strang. 




Figure 7. Peakon-antipeakon collision for 7 = 5 at time T = 
2 (left) and energy density (right) at the first time, where the 
numerical solution given by ODE45 is not positive {q = — 1.7394e — 
05). Dashdotted hne=ODE45, Stars=Exphcit Euler, Dia- 
mond=Strang. 



preserve the positivity of the particle density and conserve the invariants of our 
problem. Finally, in Figure [11] we display, with the same parameter values as 
above, the evolution in time of the energy density along the numerical solution 
given by the Strang splitting scheme. We can observe the concentration of the 
energy and then its separation in two parts, following the waves. With all these 
numerical observations, we can conclude that the proposed spatial discretisation 
is robust and qualitatively correct. The time integrators are relatively comparable 
but only the splitting schemes have the additional properties of maintaining the 
positivity of the energy density and conserve exactly the invariants of our partial 
differential equation. 
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Figure 8. Collision of smooth traveling waves: Initial data (left) 
and exact solution at time T = 11. 




Figure 9. Collision of smooth traveling waves: Evolution in time 
(left) and zoom of the evolution close to the collision. 
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Figure 10. Collision of smooth traveling waves: numerical solu- 
tions at time T = 11. Dashdotted line=ODE45, Stars=Explicit 
Euler, Diamond=Strang. 




Figure 11. Evolution of the energy density (left picture) along 
the numerical solution given by the Strang splitting and close up 
look at the blow up time (right). 
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